In [1]:
import sys
sys.path.insert(0, '/home/swl/Owncloud/Code/Edinburgh/Multiple Variables/')
# sys.path.insert(0, '/home/swk/Owncloud/Code/Edinburgh/Multiple Variables/')
In [2]:
''' Imports first etc. '''
import holoviews as hv
%load_ext holoviews.ipython
%output widgets='embed' holomap='widgets' size=200 fig='png'
import pylab as pl
import saveload as sl
import maxL as ml
from matplotlib.colors import LinearSegmentedColormap
from pylab import exp,cos,sin,pi,tan
/usr/local/lib/python2.7/dist-packages/IPython/nbformat.py:13: ShimWarning: The `IPython.nbformat` package has been deprecated. You should import from nbformat instead. "You should import from nbformat instead.", ShimWarning) /usr/local/lib/python2.7/dist-packages/IPython/nbconvert.py:13: ShimWarning: The `IPython.nbconvert` package has been deprecated. You should import from ipython_nbconvert instead. "You should import from ipython_nbconvert instead.", ShimWarning) /usr/local/lib/python2.7/dist-packages/IPython/kernel/__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead. "You should import from ipykernel or jupyter_client instead.", ShimWarning)
In [3]:
## Holoviews settings and options
# set labels
# angles = pl.linspace(-90,90,65)
# rads =
# labels =[(i,int(angles[i])) for i in pl.linspace(0,64,3)]
# dimensions
p_dim = hv.Dimension('Preferred orientation',unit='deg')
a1_dim = hv.Dimension('Stim 1',unit='deg')
a2_dim = hv.Dimension('Stim 2',unit='deg')
r_dim = hv.Dimension('Response',unit='sp/s')
L_dim = hv.Dimension('Log Likelihood')
bias_dim = hv.Dimension('Bias',unit='deg')
In [4]:
%opts Overlay [show_frame=False show_grid=False]
%opts Curve [show_grid=False]
%opts Scatter [show_grid=False show_frame=False]
%opts Raster [colorbar=True]
# yticks=labels xticks=labels] (cmap='rainbow')
Coding and decoding of multiple variables in the brain
ANC Seminar
Sander Keemink
2015-09-22
Tuning curves¶
In [5]:
reload(ml)
# function for simple firing rate
def fun(x,var=0,par = [1,1]):
''' 2 d tuning curve circular variables
A neurons tuning is simply given by multiplication of two von mises functions
Parameters
--------------
x: array
x values (preferred orientations)
var: double
orientation
par : double
Parameters. Given as [a1,k1], the inputs into von mises function.
Returns
---------------
Array of firing rates given the var and par values
'''
return ml.mises(x,var,par[0],par[1])
# set bounds for fitting
bnds = [(-pi*1.5,pi*1.5)]
# set starting conditions for fitting
var0 = [-pi/2,pi/2]
# number of values/neurons
N = 64
# number of center orientations
Nc = 16
# model vars
kc = 3
ac = 10
# basic x values for circular variables
A = pl.linspace(-pi,pi,N+1) # in radiants
Ax= A*180/pi # in degrees for plotting purposes
# centre orientation values for holomap
C = pl.linspace(-pi,pi,Nc+1)
# create empty holomaps
hmap_f = hv.HoloMap(kdims=['Centre Orientation']) # underlying firing rate
hmap_fest = hv.HoloMap(kdims=['Centre Orientation']) # underlying firing rate
hmap_r = hv.HoloMap(kdims=['Centre Orientation']) # noisy response
hmap_vl = hv.HoloMap(kdims=['Centre Orientation']) # vertical lines
# fill holomaps
for c in C: # for a
# get underlying response
f = fun(A,c,par = [ac,kc])
# get poisson response
r = pl.poisson(f)
# get decoded variables through likelihood maximization
L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par=[ac,kc],noise='poisson')
# get estimated underlying firing rate
f_est = fun(A,var=var_max,par=[ac,kc])
# add to holomap
hmap_f[c*180/pi] = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'underlying response')
hmap_fest[c*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response')
hmap_r[c*180/pi] = hv.Scatter(zip(Ax,r),vdims=[r_dim],kdims=[p_dim],label = 'measured response')
hmap_vl[c*180/pi] = hv.VLine(c*180/pi)
Coding¶
- Simple tuning curve model
In [6]:
%%opts Scatter (alpha=0) Curve (alpha=0.5) Overlay [xticks=[-180,0,180] yticks=[0,10] show_legend=False]
showfig = hmap_f*hmap_r
In [7]:
showfig
Out[7]:
In [8]:
%%opts Scatter (alpha=1) Curve (alpha=0.5) Overlay [show_legend=True]
showfig = hmap_f*hmap_r
- Measured firing rate $$ r = \texttt{Poisson}(f) $$
In [9]:
showfig
Out[9]:
- What was the encoded orientation?
Maximum Likelihood Decoding¶
Likelihood to measure a firing rate $$ p(\mathbf{r}|c) = \prod_i p(r_i|c) $$
How to decode? asymptotically optimal thing is max likelihood decoding * $$ \hat{\theta}_c = \texttt{max}_{\theta_c} p(\mathbf{r}|\theta_c) $$
*Kay (1998) - Fundamentals of Statistical Signal Processing
In [10]:
# create empty holomaps
hmap_f = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Estimated underlying firing rate
hmap_fest = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Estimated underlying firing rate
hmap_r = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Estimated underlying firing rate
hmap_L = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Likelihood
hmap_Lloc = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Likelihood
# set c
c = 0
# get underlying response
f = fun(A,c,par = [ac,kc])
curve_f = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'underlying response')
# get poisson response
r = pl.poisson(f)
scatter_r = hv.Scatter(zip(Ax,r),vdims=[r_dim],kdims=[p_dim],label = 'measured response')
# initiate empty likelihood array
L = pl.zeros(len(C))
# get Log Likelihood across estimated centre orientations
for i,c in enumerate(C): # for a
# get estimated underlying firing rate
f_est = fun(A,var=c,par=[ac,kc])
# get likelihood
L[i] = ml.findLL(f,fun,A,c,par=[ac,kc])
# add to holomap
hmap_f[c*180/pi] = curve_f
hmap_r[c*180/pi] = scatter_r
hmap_fest[c*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response')
hmap_Lloc[c*180/pi] = hv.Scatter(zip([c*180/pi],[L[i]]),vdims=[r_dim],kdims=[p_dim])(style={'s':100,'alpha':0.5})
curve_L = hv.Curve(zip(C*180/pi,L),vdims=[L_dim],kdims=[p_dim])
for i, c in enumerate(C):
hmap_L[c*180/pi] = curve_L
In [11]:
%%opts Scatter (alpha=1) Curve (alpha=0.5) Overlay [xticks=[-180,0,180] yticks=3 show_legend=True]
showfig = hmap_f*hmap_r*hmap_fest+hmap_L*hmap_Lloc*hmap_Lloc
In [12]:
showfig
Out[12]:
In [13]:
%%opts Scatter (alpha=1) Curve (alpha=0.5) Overlay [xticks=[-180,0,180] yticks=[0,10] show_legend=True]
showfig = hmap_f*hmap_fest*hmap_r
In [14]:
showfig
Out[14]:
- 'Optimal' variance (cramer-rao bound)
- Non biased
- I will be focussing on the bias
But neurons rarely respond to only single stimulus¶
- overlapping stimuli
- other variables/contexts (color, motion, speed, etc.)
Fig: Schwartz, Sejnowski & Dayan (2009)
Adjust coding model¶
$$f(\theta_1,\theta_2) = a\exp(k_1\cos[\phi-\theta_1])\{1-b\exp(k_2\cos[\phi-\theta_2])\} $$In [15]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
''' 2 d tuning curve circular variables
A neurons tuning is simply given by multiplication of two von mises functions
Parameters
--------------
x: array
x values (preferred orientations)
var: array
value of the two inputs
par : array
Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
Returns
---------------
Array of firing rates given the var and par values
'''
return ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(x,var[1],par[2],par[3]))
# set parameters
par = [10,3,0.5,2]
par1 = [10,3,0,1] # effect of second stimulus turned off
par2 = [1,0,0.5,1] # effect of first stimulus turned off
# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))
# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]
# number of neurons
N = 32
# number of contexts
nA = 16
# iterations
iterations = 500
# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))
# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)
# stim 1
a1 = 0
# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)
# predefine holomaps
hmap_scatters = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1 = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2 = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f = hv.HoloMap(kdims=['context']) # population response
hmap_g = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h = hv.HoloMap(kdims=['context']) # stim 2 modulation
# loop over stim 2's
for i,a2 in enumerate(A2):
print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
# set variables
var=[a1,a2]
# get underlying firing rate
f = fun(A,var=var,par=par)
hmap_f[a2*180/pi] = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'response')
# get drive and modulation
g = fun(A,var=var,par=par1)
hmap_g[a2*180/pi] = hv.Curve(zip(Ax,g),vdims=[r_dim],kdims=[p_dim],label = 'Stim 1 drive')
h = fun(A,var=var,par=par2)
hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
# do noise realizations
print 'Realization (max: ' + str(iterations) + '):'
for it in range(iterations):
print '\r'+str(it),
# get noisy firing rate
r = pl.poisson(f)
# do the maximization
L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
# store variables
a1_est[i,it] = var_max[0]
a2_est[i,it] = var_max[1]
# make sure everything is within -pi,pi
a1_est[i,a1_est[i,:]<-pi]+=2*pi
a1_est[i,a1_est[i,:]>pi] -=2*pi
a2_est[i,a2_est[i,:]<-pi]+=2*pi
a2_est[i,a2_est[i,:]>pi] -=2*pi
# get bias
a1_bias[i] = ml.circstats(a1_est[i,:])[0]
a2_meanest = ml.circstats(a2_est[i,:])[0]
# compare estimate for a2 with actual a2 to make circular
a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
a2_bias[a2_bias<-pi]+=2*pi
a2_bias[a2_bias>pi] -=2*pi
# fill in relevant holomaps
hmap_scatters[a2*180/pi] = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
# fill in illusion curves
for i,a2 in enumerate(A2):
hmap_illusion1[a2*180/pi] = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
hmap_illusion2[a2*180/pi] = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16 Realization (max: 500): 499 Doing stimulus 1 of 16 Realization (max: 500): 499 Doing stimulus 2 of 16 Realization (max: 500): 499 Doing stimulus 3 of 16 Realization (max: 500): 499 Doing stimulus 4 of 16 Realization (max: 500): 499 Doing stimulus 5 of 16 Realization (max: 500): 499 Doing stimulus 6 of 16 Realization (max: 500): 499 Doing stimulus 7 of 16 Realization (max: 500): 499 Doing stimulus 8 of 16 Realization (max: 500): 499 Doing stimulus 9 of 16 Realization (max: 500): 499 Doing stimulus 10 of 16 Realization (max: 500): 499 Doing stimulus 11 of 16 Realization (max: 500): 499 Doing stimulus 12 of 16 Realization (max: 500): 499 Doing stimulus 13 of 16 Realization (max: 500): 499 Doing stimulus 14 of 16 Realization (max: 500): 499 Doing stimulus 15 of 16 Realization (max: 500): 499
In [16]:
showfig = hmap_f*hmap_g*hmap_h
In [17]:
showfig
Out[17]:
In [18]:
%%opts Curve {+axiswise} Scatter {+axiswise}
showfig = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g*hmap_h
In [19]:
showfig.cols(2)
Out[19]:
In [20]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g*hmap_h
showfig += hmap_scatters
In [21]:
showfig.cols(2)
Out[21]:
In [22]:
labels =[(i,int(A[i]*180/pi)) for i in pl.linspace(0,len(A)-1,3)]
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:1: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future if __name__ == '__main__':
In [23]:
%%opts Raster [colorbar=True yticks=labels xticks=labels] (cmap='rainbow')
%%opts Scatter (alpha=0.25 color='k')
# neruongs
nA = 32
A = pl.arange(-pi,pi,2*pi/nA)
# sampling
nS = 32
As = pl.arange(-pi,pi,2*pi/nS)
# set variables
varid = 11
var = [0,A2[varid]] #a1,a2
# get underlying firing rate
f = fun(A,var=var,par=par)
curve_f = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'underlying response')
# predefine holomaps
hmap_fest1= hv.HoloMap(kdims=['stim2']) # underlying population response
hmap_fest2= hv.HoloMap(kdims=['stim2']) # underlying population response
hmap_f = hv.HoloMap(kdims=['stim2']) # estimated population response
hmap_pos1 = hv.HoloMap(kdims=['stim2']) # position of stimulus
hmap_pos2 = hv.HoloMap(kdims=['stim2']) # position of stimulus
hmap_L = hv.HoloMap(kdims=['stim2']) # likelihood map
# get likelihood across estimations
L = pl.zeros((len(A),len(A)))
for i,a1 in enumerate(A):
for j,a2 in enumerate(A):
# get likelihood
L[i,j] = ml.findLL(f,fun,A,[a1,a2],par)
a1a = 0
a1b = -pi/16
for i,a2 in enumerate(As):
hmap_f[a2*180/pi] = curve_f
# get underlying firing rate situation 1
f_est = fun(A,var=[a1a,a2],par=par)
hmap_fest1[a2*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response 1')
# situation 2
f_est = fun(A,var=[a1b,a2],par=par)
hmap_fest2[a2*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response 2')
# get normalized version of relevant scatter
var1_norm = len(A)*(a1_est[varid,:]+pi)/(2*pi)
var2_norm = len(A)*(a2_est[varid,:]+pi)/(2*pi)
scatter = hv.Scatter(zip(var2_norm,var1_norm))
showfig = hmap_f*hmap_fest1*hmap_fest2
Ambiguity in encoding model¶
In [24]:
showfig
Out[24]:
In [25]:
showfig = hv.Raster(L,kdims=[a2_dim,a1_dim])*scatter+hv.Empty()
In [26]:
showfig
Out[26]:
Stimulus-1-dependence¶
- from
- to
In [27]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
''' 2 d tuning curve circular variables
A neurons tuning is simply given by multiplication of two von mises functions
Parameters
--------------
x: array
x values (preferred orientations)
var: array
value of the two inputs
par : array
Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
Returns
---------------
Array of firing rates given the var and par values
'''
return ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(var[0],var[1],par[2],par[3]))
# set parameters
par = [10,3,0.5,2]
par1 = [10,3,0,1] # effect of second stimulus turned off
par2 = [1,0,0.5,1] # effect of first stimulus turned off
# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))
# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]
# number of neurons
N = 32
# number of contexts
nA = 16
# iterations
iterations = 500
# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))
# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)
# stim 1
a1 = 0
# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)
# predefine holomaps
hmap_scatters = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1 = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2 = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f = hv.HoloMap(kdims=['context']) # population response
hmap_g = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h = hv.HoloMap(kdims=['context']) # stim 2 modulation
# loop over stim 2's
for i,a2 in enumerate(A2):
print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
# set variables
var=[a1,a2]
# get underlying firing rate
f = fun(A,var=var,par=par)
hmap_f[a2*180/pi] = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'response')
# get drive and modulation
g = fun(A,var=var,par=par1)
hmap_g[a2*180/pi] = hv.Curve(zip(Ax,g),vdims=[r_dim],kdims=[p_dim],label = 'Stim 1 drive')
h = fun(A,var=var,par=par2)
hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
# do iterations
print 'iteration (max: ' + str(iterations) + '):'
for it in range(iterations):
print '\r'+str(it),
# get noisy firing rate
r = pl.poisson(f)
# do the maximization
L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
# store variables
a1_est[i,it] = var_max[0]
a2_est[i,it] = var_max[1]
# by convenction set sign properly (since can't be told by decoder)
# which
if pl.sign(a2_est[i,it]) != pl.sign(a2):
a2_est[i,it] *= -1
# make sure everything is within -pi,pi
a1_est[i,a1_est[i,:]<-pi]+=2*pi
a1_est[i,a1_est[i,:]>pi] -=2*pi
a2_est[i,a2_est[i,:]<-pi]+=2*pi
a2_est[i,a2_est[i,:]>pi] -=2*pi
# get bias
a1_bias[i] = ml.circstats(a1_est[i,:])[0]
a2_meanest = ml.circstats(a2_est[i,:])[0]
# compare estimate for a2 with actual a2 to make circular
a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
a2_bias[a2_bias<-pi]+=2*pi
a2_bias[a2_bias>pi] -=2*pi
# fill in relevant holomaps
hmap_scatters[a2*180/pi] = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
# fill in illusion curves
for i,a2 in enumerate(A2):
hmap_illusion1[a2*180/pi] = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
hmap_illusion2[a2*180/pi] = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16 iteration (max: 500): 499 Doing stimulus 1 of 16 iteration (max: 500): 499 Doing stimulus 2 of 16 iteration (max: 500): 499 Doing stimulus 3 of 16 iteration (max: 500): 499 Doing stimulus 4 of 16 iteration (max: 500): 499 Doing stimulus 5 of 16 iteration (max: 500): 499 Doing stimulus 6 of 16 iteration (max: 500): 499 Doing stimulus 7 of 16 iteration (max: 500): 499 Doing stimulus 8 of 16 iteration (max: 500): 499 Doing stimulus 9 of 16 iteration (max: 500): 499 Doing stimulus 10 of 16 iteration (max: 500): 499 Doing stimulus 11 of 16 iteration (max: 500): 499 Doing stimulus 12 of 16 iteration (max: 500): 499 Doing stimulus 13 of 16 iteration (max: 500): 499 Doing stimulus 14 of 16 iteration (max: 500): 499 Doing stimulus 15 of 16 iteration (max: 500): 499
In [28]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g*hmap_h
showfig += hmap_scatters
In [29]:
showfig.cols(2)
Out[29]:
Multiple populations¶
$$f_a(\theta_1,\theta_2) = a_a\exp(k_{a1}\cos[\phi-\theta_1])\{1-b\exp(k_{a2}\cos[\phi-\theta_2])\} $$$$f_b(\theta_1,\theta_2) = a_b\exp(k_{b1}\cos[\phi-\theta_2])\{1-b\exp(k_{b2}\cos[\phi-\theta_1])\} $$In [30]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
''' 2 d tuning curve circular variables
A neurons tuning is simply given by multiplication of two von mises functions
Parameters
--------------
x: array
x values (preferred orientations)
var: array
value of the two inputs
par : array
Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
Returns
---------------
Array of firing rates given the var and par values
'''
f = pl.zeros(2*len(x))
f[:len(x)] = ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(x,var[1],par[2],par[3]))
f[len(x):] = ml.mises(x,var[1],par[0],par[1])*(1-ml.mises(x,var[0],par[2],par[3]))
return f
# set parameters
par = [10,3,0.5,2,10,3,0.5,2]
# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))
# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]
# number of neurons
N = 32
# basic x values for circular variables
A = pl.linspace(-pi,pi,N+1) # in radiants
Ax = pl.linspace(-pi,3*pi,2*N+1)
Ax= Ax*180/pi # in degrees for plotting purposes
# number of contexts
nA = 16
# iterations
iterations = 500
# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))
# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)
# stim 1
a1 = 0
# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)
# predefine holomaps
hmap_scatters = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1 = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2 = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f = hv.HoloMap(kdims=['context']) # population response
hmap_g = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h = hv.HoloMap(kdims=['context']) # stim 2 modulation
# loop over stim 2's
for i,a2 in enumerate(A2):
print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
# set variables
var=[a1,a2]
# get underlying firing rate
f = fun(A,var=var,par=par)
hmap_f[a2*180/pi] = hv.Curve(zip(Ax[:N],f[:N]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 1')
# get second pop
# g = fun(A,var=var,par=par1)
hmap_g[a2*180/pi] = hv.Curve(zip(Ax[N:],f[N:]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 2')
# get modulation
h = fun(A,var=var,par=par2)
hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
# do iterations
print 'iteration (max: ' + str(iterations) + '):'
for it in range(iterations):
print '\r'+str(it),
# get noisy firing rate
r = pl.poisson(f)
# do the maximization
L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
# store variables
a1_est[i,it] = var_max[0]
a2_est[i,it] = var_max[1]
# by convenction set sign properly (since can't be told by decoder)
# which
# if pl.sign(a2_est[i,it]) != pl.sign(a2):
# a2_est[i,it] *= -1
# make sure everything is within -pi,pi
a1_est[i,a1_est[i,:]<-pi]+=2*pi
a1_est[i,a1_est[i,:]>pi] -=2*pi
a2_est[i,a2_est[i,:]<-pi]+=2*pi
a2_est[i,a2_est[i,:]>pi] -=2*pi
# get bias
a1_bias[i] = ml.circstats(a1_est[i,:])[0]
a2_meanest = ml.circstats(a2_est[i,:])[0]
# compare estimate for a2 with actual a2 to make circular
a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
a2_bias[a2_bias<-pi]+=2*pi
a2_bias[a2_bias>pi] -=2*pi
# fill in relevant holomaps
hmap_scatters[a2*180/pi] = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
# fill in illusion curves
for i,a2 in enumerate(A2):
hmap_illusion1[a2*180/pi] = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
hmap_illusion2[a2*180/pi] = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16 iteration (max: 500): 499 Doing stimulus 1 of 16 iteration (max: 500): 499 Doing stimulus 2 of 16 iteration (max: 500): 499 Doing stimulus 3 of 16 iteration (max: 500): 499 Doing stimulus 4 of 16 iteration (max: 500): 499 Doing stimulus 5 of 16 iteration (max: 500): 499 Doing stimulus 6 of 16 iteration (max: 500): 499 Doing stimulus 7 of 16 iteration (max: 500): 499 Doing stimulus 8 of 16 iteration (max: 500): 499 Doing stimulus 9 of 16 iteration (max: 500): 499 Doing stimulus 10 of 16 iteration (max: 500): 499 Doing stimulus 11 of 16 iteration (max: 500): 499 Doing stimulus 12 of 16 iteration (max: 500): 499 Doing stimulus 13 of 16 iteration (max: 500): 499 Doing stimulus 14 of 16 iteration (max: 500): 499 Doing stimulus 15 of 16 iteration (max: 500): 499
In [31]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g#*hmap_h
showfig += hmap_scatters
In [32]:
showfig.cols(2)
Out[32]:
Asymetric populations¶
In [33]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
''' 2 d tuning curve circular variables
A neurons tuning is simply given by multiplication of two von mises functions
Parameters
--------------
x: array
x values (preferred orientations)
var: array
value of the two inputs
par : array
Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
Returns
---------------
Array of firing rates given the var and par values
'''
f = pl.zeros(2*len(x))
f[:len(x)] = ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(x,var[1],par[2],par[3]))
f[len(x):] = ml.mises(x,var[1],par[4],par[5])*(1-ml.mises(x,var[0],par[6],par[7]))
return f
# set parameters
par = [10,3,0.5,2,8,1.5,0.7,1]
# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))
# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]
# number of neurons
N = 32
# basic x values for circular variables
A = pl.linspace(-pi,pi,N+1) # in radiants
Ax = pl.linspace(-pi,3*pi,2*N+1)
Ax= Ax*180/pi # in degrees for plotting purposes
# number of contexts
nA = 16
# iterations
iterations = 500
# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))
# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)
# stim 1
a1 = 0
# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)
# predefine holomaps
hmap_scatters = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1 = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2 = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f = hv.HoloMap(kdims=['context']) # population response
hmap_g = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h = hv.HoloMap(kdims=['context']) # stim 2 modulation
# loop over stim 2's
for i,a2 in enumerate(A2):
print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
# set variables
var=[a1,a2]
# get underlying firing rate
f = fun(A,var=var,par=par)
hmap_f[a2*180/pi] = hv.Curve(zip(Ax[:N],f[:N]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 1')
# get second pop
# g = fun(A,var=var,par=par1)
hmap_g[a2*180/pi] = hv.Curve(zip(Ax[N:],f[N:]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 2')
# get modulation
# h = fun(A,var=var,par=par2)
# hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
# do iterations
print 'iteration (max: ' + str(iterations) + '):'
for it in range(iterations):
print '\r'+str(it),
# get noisy firing rate
r = pl.poisson(f)
# do the maximization
L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
# store variables
a1_est[i,it] = var_max[0]
a2_est[i,it] = var_max[1]
# by convenction set sign properly (since can't be told by decoder)
# which
# if pl.sign(a2_est[i,it]) != pl.sign(a2):
# a2_est[i,it] *= -1
# make sure everything is within -pi,pi
a1_est[i,a1_est[i,:]<-pi]+=2*pi
a1_est[i,a1_est[i,:]>pi] -=2*pi
a2_est[i,a2_est[i,:]<-pi]+=2*pi
a2_est[i,a2_est[i,:]>pi] -=2*pi
# get bias
a1_bias[i] = ml.circstats(a1_est[i,:])[0]
a2_meanest = ml.circstats(a2_est[i,:])[0]
# compare estimate for a2 with actual a2 to make circular
a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
a2_bias[a2_bias<-pi]+=2*pi
a2_bias[a2_bias>pi] -=2*pi
# fill in relevant holomaps
hmap_scatters[a2*180/pi] = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
# fill in illusion curves
for i,a2 in enumerate(A2):
hmap_illusion1[a2*180/pi] = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
hmap_illusion2[a2*180/pi] = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16 iteration (max: 500): 499 Doing stimulus 1 of 16 iteration (max: 500): 499 Doing stimulus 2 of 16 iteration (max: 500): 499 Doing stimulus 3 of 16 iteration (max: 500): 499 Doing stimulus 4 of 16 iteration (max: 500): 499 Doing stimulus 5 of 16 iteration (max: 500): 499 Doing stimulus 6 of 16 iteration (max: 500): 499 Doing stimulus 7 of 16 iteration (max: 500): 499 Doing stimulus 8 of 16 iteration (max: 500): 499 Doing stimulus 9 of 16 iteration (max: 500): 499 Doing stimulus 10 of 16 iteration (max: 500): 499 Doing stimulus 11 of 16 iteration (max: 500): 499 Doing stimulus 12 of 16 iteration (max: 500): 499 Doing stimulus 13 of 16 iteration (max: 500): 499 Doing stimulus 14 of 16 iteration (max: 500): 499 Doing stimulus 15 of 16 iteration (max: 500): 499
In [34]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g#*hmap_h
showfig += hmap_scatters
In [35]:
showfig.cols(2)
Out[35]:
Conclusion/discussion¶
- Multiple variables make decoding/coding a lot harder due to ambiguities
- Biases and decoding distributions can be very non-trivial
- Could make predictions about bias distributions in humans
Thanks to...¶
- Mark van Rossum, supervisor
- Holoviews
Notes to self¶
- multiple variables are okay, if other variables are for same preferred variable across population